#-----------------------------------------------------------------------------#
#### 0. IMPORT ####
#-----------------------------------------------------------------------------#

gc()
memory.limit(size=64000)

library(data.table)
library("lfe")
library(texreg)
library(readstata13)
# library(plm)
library(Hmisc)
library(haven)
library(plyr)

# Functions
simplelag <- function (x,by=NULL,mylag=1,outside=NA) {
  myend <- length(x)- mylag
  if(!is.null(by))  {
    lby <- c(replicate(mylag,""),as.character(by[1:myend]))
    y0 <- c(replicate(mylag,outside),x[1:myend])
    y <- ifelse(as.character(by)==lby,y0,outside)
  }
  else {
    y <- c(replicate(mylag,outside),x[1:myend])
  }
}

retain <- function(x,event,outside=NA)
{
  indices <- c(1,which(event==TRUE),length(x)+1)
  values <- c(outside,x[event %in% TRUE])
  y <- rep(values,diff(indices))
}

obtain <- function(x,event,outside=NA)
{
  indices <- c(0,which(event==TRUE),length(x))
  values <- c(x[event %in% TRUE],outside)
  y <- rep(values,diff(indices))
}


wtd.summary <- function (x,w=NA) {
  name <- deparse(substitute(x))
  nb_obs <- length(x)
  sum_wgt <- sum(w[is.na(x) ==F])
  nb_missing <- sum(is.na(x))
  sum_wgt_missing <- sum(w[is.na(x) ==T])
  wtd_mean <- weighted.mean(x,w=w,na.rm=T)
  wtd_sd <- wtd.var(x,w=w,na.rm=T)**0.5
  min <- min(x,na.rm=T)
  q1 <- wtd.quantile(x,weights=w,probs=0.25)
  q2 <- wtd.quantile(x,weights=w,probs=0.50)
  q3 <- wtd.quantile(x,weights=w,probs=0.75)
  max <- max(x,na.rm=T)
  structure(data.frame(name,nb_obs,sum_wgt,nb_missing,sum_wgt_missing,wtd_mean,wtd_sd,min,q1,q2,q3,max))
}
setwd("C:\\Users\\Public\\Documents\\Olivier\\Results\\Establishment_Composition")



#DATA
ee <- readRDS("outsourcing.rds")
str(ee)

ee[,c("sf9910","f9910_w","f9910xf9910",
      "f9099","f9910","f9010","sf9010","lndnbwkrs_neg",
      "lcount","lnbwkrs","one"):=NULL] #Deleting variables

ee <- ee[ee$f9010_w>0 & is.na(ee$f9010_w)==F,]

ee$routs_lsk <- ifelse(is.na(ee$depsize)==T | ee$depsize==0,0,ee$outs_lsk/ee$depsize)
ee$routs_arr_lsk <- ifelse(is.na(ee$deparrsize)==T | ee$deparrsize==0,
                           ifelse(is.na(ee$arrdepsize)==T | ee$arrdepsize==0,0,ee$arr_lsk/ee$arrdepsize),
                           ee$outs_lsk/ee$deparrsize)
ee$yrouts_lsk <- ave(ee$routs_lsk*(ee$youts_lsk==10000 | ee$year<ee$y2outs_lsk),ee$est,FUN=max)
ee$yrouts_arr_lsk <- ave(ee$routs_arr_lsk*(pmin(ee$youts_lsk,ee$yarr_lsk)==10000 | ee$year<pmin(ee$y2outs_lsk,ee$y2arr_lsk)),ee$est,FUN=max)

  # ee <- ee[order(ee$est,ee$year),]
  # rownames(ee) <- NULL
  # ee$nouts_lsk <- ave((ee$outs_lsk>0),ee$est,FUN=sum)
  # ee$cumouts_lsk <- ave((ee$outs_lsk>0),ee$est,FUN=cumsum)
  # There's 2664 1st events, and only 44 2nd events and 11 3rd events. Hence this justify focusing on 1st events only
  # ee[,c("cumouts_lsk"):=NULL] #Deleting variables
  
  # ee$y2outs_lsk <- ave(ifelse(ee$year==ee$youts_lsk,10000,ifelse((ee$outs_lsk>0)*ee$year>0,(ee$outs_lsk>0)*ee$year,10000)),ee$est,FUN=min)
  
str(ee)
gc()

table(ee$year)
table(ee$youts_lsk)
str(ee)
gc()

#-----------------------------------------------------------------------------#
#### 1.OUTSOURCING DATA ####
#-----------------------------------------------------------------------------#


fff <- NULL
for (y in c(2002:2017)) {
ff <- ee[ee$year>y-5 & ee$year<y+5
         & (ee$youts_lsk %in% c(y,10000) | ee$year<ee$youts_lsk)
         & ee$year<ee$y2outs_lsk,
                 c("f9010xf9010","f9010_w","year","youts_lsk","yrouts_lsk","est","count","lnnbwkrs_cumneg")]
ff$eyear <- y
fff <- rbind(fff,ff)
}
rm(ff)
gc()
nrow(fff)

set.seed(123)

rrr <- unique(fff[,c("eyear","est")])
rrr$rand <- runif(nrow(rrr))
fff <- merge(fff,rrr,by=c("eyear","est"))
str(fff)
rm(rrr)
fff <- fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,]
fff$firm <- substr(fff$est,1,9)

#-----------------------------------------------------------------------------#
#### 2.OUTSOURCING ESTIMATES ####
#-----------------------------------------------------------------------------#

# EVENT
gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
            + I((year - youts_lsk) %in% -4 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% -3 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% -2 & youts_lsk==eyear)
            # + I((year - youts_lsk) %in% -1 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% 0 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% 1 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% 2 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% 3 & youts_lsk==eyear)
            + I((year - youts_lsk) %in% 4 & youts_lsk==eyear)
            |paste(eyear,est)|0|firm, 
            weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
            data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_lsk_stacked.html")


# EVENT + SIZE CTRL
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I((year - youts_lsk) %in% -4 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% -3 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% -2 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% -1 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 0 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 1 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 2 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 3 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 4 & youts_lsk==eyear)
           +log(count)
           +lnnbwkrs_cumneg
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_lsk_stacked_size.html")
str(fff)

# EVENT MAGNITUDE
gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I(yrouts_lsk*((year - youts_lsk) %in%  -4 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% -3 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% -2 & youts_lsk==eyear))
           # + I(yrouts_lsk*((year - youts_lsk) %in% -1 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 0 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 1 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 2 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 3 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk)  %in% 4 & youts_lsk==eyear))
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_rlsk_stacked.html")


# EVENT MAGNITUDE + SIZE CTRL
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I(yrouts_lsk*((year - youts_lsk) %in%  -4 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% -3 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% -2 & youts_lsk==eyear))
           # + I(yrouts_lsk*((year - youts_lsk) %in% -1 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 0 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 1 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 2 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 3 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk)  %in% 4 & youts_lsk==eyear))
           +log(count)
           +lnnbwkrs_cumneg
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_rlsk_stacked_size.html")
str(fff)

# Trend
# EVENT

gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ year
           # + I((year - youts_lsk) %in% -4 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% -3 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% -2 & youts_lsk==eyear)
           # # + I((year - youts_lsk) %in% -1 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% 0 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% 1 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% 2 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% 3 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% 4 & youts_lsk==eyear)
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01),digits=5)
htmlreg(list(dd),stars=c(0.1,0.05,0.01),digits=5,file="twfe_stacked_trend.html")

rm(dd)
dd <- felm(I(100*f9010xf9010) ~ year
           + I((year - youts_lsk) %in% -4 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% -3 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% -2 & youts_lsk==eyear)
           # + I((year - youts_lsk) %in% -1 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 0 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 1 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 2 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 3 & youts_lsk==eyear)
           + I((year - youts_lsk) %in% 4 & youts_lsk==eyear)
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01),digits=5)
htmlreg(list(dd),stars=c(0.1,0.05,0.01),digits=5,file="twfe_lsk_stacked_trend1.html")


gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ year+
           + I(yrouts_lsk*((year - youts_lsk) %in%  -4 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% -3 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% -2 & youts_lsk==eyear))
           # + I(yrouts_lsk*((year - youts_lsk) %in% -1 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 0 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 1 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 2 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk) %in% 3 & youts_lsk==eyear))
           + I(yrouts_lsk*((year - youts_lsk)  %in% 4 & youts_lsk==eyear))
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01),digits=5)
htmlreg(list(dd),stars=c(0.1,0.05,0.01),digits=5,file="twfe_rlsk_stacked_trend.html")



#-----------------------------------------------------------------------------#
#### 3.DEP+ARR OUT DATA ####
#-----------------------------------------------------------------------------#
rm(fff)
ee$youts_arr_lsk <- pmin(ee$youts_lsk,ee$yarr_lsk)
ee$y2outs_arr_lsk <- pmin(ee$y2outs_lsk,ee$y2arr_lsk)

gc()
fff <- NULL
for (y in c(2002:2019)) {
  ff <- ee[ee$year>y-5 & ee$year<y+5
           & (ee$youts_arr_lsk %in% c(y,10000) | ee$year<ee$youts_arr_lsk)
           & ee$year<ee$y2outs_arr_lsk,
           c("f9010xf9010","f9010_w","year","youts_arr_lsk","yrouts_arr_lsk","est","count","lnnbwkrs_cumneg")]
  ff$eyear <- y
  fff <- rbind(fff,ff)
}
rm(ff)
gc()
nrow(fff)

set.seed(123)

rrr <- unique(fff[,c("eyear","est")])
rrr$rand <- runif(nrow(rrr))
fff <- merge(fff,rrr,by=c("eyear","est"))
fff <- fff[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9,]
fff$firm <- substr(fff$est,1,9)
str(fff)
rm(rrr)

#-----------------------------------------------------------------------------#
#### 4.DEP+ARR OUT. EST. ####
#-----------------------------------------------------------------------------#

gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I((year - youts_arr_lsk) %in% -4 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% -3 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% -2 & youts_arr_lsk==eyear)
           # + I((year - youts_arr_lsk) %in% -1 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 0 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 1 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 2 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 3 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 4 & youts_arr_lsk==eyear)
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_lsk_deparr_stacked.html")


gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I((year - youts_arr_lsk) %in% -4 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% -3 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% -2 & youts_arr_lsk==eyear)
           # + I((year - youts_arr_lsk) %in% -1 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 0 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 1 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 2 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 3 & youts_arr_lsk==eyear)
           + I((year - youts_arr_lsk) %in% 4 & youts_arr_lsk==eyear)
           + log(count)
           + lnnbwkrs_cumneg
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_lsk_deparr_stacked_size.html")

gc()
rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -4 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -3 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -2 & youts_arr_lsk==eyear))
           # + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -1 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 0 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 1 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 2 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 3 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 4 & youts_arr_lsk==eyear))
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_rlsk_deparr_stacked.html")


rm(dd)
gc()
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -4 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -3 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -2 & youts_arr_lsk==eyear))
           # + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% -1 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 0 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 1 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 2 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 3 & youts_arr_lsk==eyear))
           + I(yrouts_arr_lsk*((year - youts_arr_lsk) %in% 4 & youts_arr_lsk==eyear))
           +log(count)
           +lnnbwkrs_cumneg
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9],
           data=fff[fff$eyear==fff$youts_arr_lsk | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(list(dd),stars=c(0.1,0.05,0.01),file="twfe_rlsk_deparr_stacked_size.html")





#-----------------------------------------------------------------------------#
#### 3.DEP+ARR SUB DATA ####
#-----------------------------------------------------------------------------#
rm(fff)
rm(dd)

ee$ydep_sub <- ave(ifelse((ee$depflow_sub>0)*ee$year>0,(ee$depflow_sub>0)*ee$year,10000),ee$est,FUN=function(x) min(x,na.rm=T))
ee$yarr_sub <- ave(ifelse((ee$arrflow_sub>0)*ee$year>0,(ee$arrflow_sub>0)*ee$year,10000),ee$est,FUN=function(x) min(x,na.rm=T))
ee$y2arr_sub <- ave(ifelse(ee$year==ee$yarr_sub,10000,ifelse((ee$arrflow_sub>0)*ee$year>0,(ee$arrflow_sub>0)*ee$year,10000)),ee$est,FUN=min)
ee$y2dep_sub <- ave(ifelse(ee$year==ee$ydep_sub,10000,ifelse((ee$depflow_sub>0)*ee$year>0,(ee$depflow_sub>0)*ee$year,10000)),ee$est,FUN=min)
ee$rdep_sub <- ifelse(is.na(ee$depsize)==T | ee$depsize==0,0,ee$depflow_sub/ee$depsize)
ee$rdep_arr_sub <- ifelse(is.na(ee$deparrsize)==T | ee$deparrsize==0,
                          ifelse(is.na(ee$arrdepsize)==T | ee$arrdepsize==0,0,ee$arrflow_sub/ee$arrdepsize),
                          ee$depflow_sub/ee$deparrsize)
ee$yrdep_sub <- ave(ee$rdep_sub*(ee$ydep_sub==10000 | ee$year<ee$y2dep_sub),ee$est,FUN=max)
ee$yrdep_arr_sub <- ave(ee$rdep_arr_sub*(pmin(ee$ydep_sub,ee$yarr_sub)==10000 | ee$year<pmin(ee$y2dep_sub,ee$y2arr_sub)),ee$est,FUN=max)

ee$ydep_arr_sub <- pmin(ee$ydep_sub,ee$yarr_sub)
ee$y2dep_arr_sub <- pmin(ee$y2dep_sub,ee$y2arr_sub)



gc()
fff <- NULL
for (y in c(2002:2019)) {
  ff <- ee[ee$year>y-5 & ee$year<y+5
           & (ee$ydep_arr_sub %in% c(y,10000) | ee$year<ee$ydep_arr_sub)
           & ee$year<ee$y2dep_arr_sub,
           c("f9010xf9010","f9010_w","year","ydep_arr_sub","yrdep_arr_sub","est","count","lnnbwkrs_cumneg")]
  ff$eyear <- y
  fff <- rbind(fff,ff)
}
rm(ff)
gc()
nrow(fff)

set.seed(123)

rrr <- unique(fff[,c("eyear","est")])
rrr$rand <- runif(nrow(rrr))
fff <- merge(fff,rrr,by=c("eyear","est"))
fff <- fff[fff$eyear==fff$ydep_arr_sub | fff$rand>0.9,]
fff$firm <- substr(fff$est,1,9)

str(fff)
rm(rrr)

rm(dd)
gc()
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I((((year - ydep_arr_sub) %in% c(-4)) & ydep_arr_sub==eyear))
           + I(( ((year - ydep_arr_sub) %in% -3)  & ydep_arr_sub==eyear))
           + I(( ((year - ydep_arr_sub) %in% -2)  & ydep_arr_sub==eyear))
           # + I(( ((year - ydep_sub) %in% -1) ))
           + I(( ((year - ydep_arr_sub) %in% 0)  & ydep_arr_sub==eyear))
           + I(( ((year - ydep_arr_sub) %in% 1) & ydep_arr_sub==eyear))
           + I(( ((year - ydep_arr_sub) %in% 2) & ydep_arr_sub==eyear))
           + I(( ((year - ydep_arr_sub) %in% 3)  & ydep_arr_sub==eyear))
           + I(( ((year - ydep_arr_sub)  %in% 4) & ydep_arr_sub==eyear))
           |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$ydep_arr_sub | fff$rand>0.9],
           data=fff[fff$eyear==fff$ydep_arr_sub | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(dd,stars=c(0.1,0.05,0.01),file="twfe_sub_deparr_stacked.html")

rm(dd)
dd <- felm(I(100*f9010xf9010) ~ factor(paste(eyear,year))
           + I(yrdep_arr_sub*(((year - ydep_arr_sub) %in% c(-4)) & ydep_arr_sub==eyear))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub) %in% -3)  & ydep_arr_sub==eyear))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub) %in% -2)  & ydep_arr_sub==eyear))
           # + I(yrdep_arr_sub*( ((year - ydep_sub) %in% -1) ))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub) %in% 0)  & ydep_arr_sub==eyear))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub) %in% 1) & ydep_arr_sub==eyear))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub) %in% 2) & ydep_arr_sub==eyear))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub) %in% 3)  & ydep_arr_sub==eyear))
           + I(yrdep_arr_sub*( ((year - ydep_arr_sub)  %in% 4) & ydep_arr_sub==eyear))
              |paste(eyear,est)|0|firm, 
           weights=fff$f9010_w[fff$eyear==fff$ydep_arr_sub | fff$rand>0.9],
           data=fff[fff$eyear==fff$ydep_arr_sub | fff$rand>0.9,])
screenreg(list(dd),stars=c(0.1,0.05,0.01))
htmlreg(dd,stars=c(0.1,0.05,0.01),file="twfe_rsub_deparr_stacked.html")

gc()
